library(parallel)
library(doMC)
library(DHARMa)
library(broom.mixed)
library(MuMIn)
library(mgcv)
library(bbmle)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2)
library(tidyr)
library(purrr)
library(plyr)
library(dplyr)

1 - Variables

table 1: Variáveis utilizadas na descrição estatística

Variables
code description class range
n_nRef number of SADs not refuted by the KS bootstrap test with critical p-value 0.05 continuous (0 ; 100)
diffS_mean diffS = (S_MN - S_obs) / S_obs, diffS_mean = mean(diffS); S_MN = species richness estimated by MN continuous (-0.977 ; 4.78)
U_med average of 10 replicates of the speciation rate estimated by MNEE continuous (8.860e-5 ; 4.221e-2)
p proportion of tree cover continuous (0.013 ; 0.961)
MN Neutral Model (EE = spatial explicit; EI = spatial implict) category 2 levels
d mean dispersal distance (meters) continuous (0.456 ; 19.183)
k proportion of propagules until the nearest neighborhood seq(0.99 : 0.05) category 20 levels
d_Lplot mean dispersal distance / side of the sample area continuous (0.02 ; 0.192)
S_obs observed species richness integer (26 ; 458)
Ntotal number of individuals in the sample area integer (649 ; 12105)
SiteCode sample area code category 103 levels

Figure 1 Possible Predictor Variables

2 - Number of unrefuted SADs

Figura 2.1 número de SADs não refutadas ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.

Figura 2.2 número de SADs não refutadas ~ (d * MN|Site ~ p). Site está ordenado pelo valor de p.

2.1 GLMM binomial

Modelo cheio 1:

linear term

  1. ~ (p + p^2) * (d + d^2) * MN; if : var_dispersao = contiguous
  2. ~ ( p + p^2) * k * MN; if : var_dispersao = category

random term

  1. 1|SiteCode
  2. MN|SiteCode
  3. (var_dispersao + var_dispersao^2)*MN|SiteCode, if : var = dispersao : contiguous
# dados
df_md <- df_resultados %>% 
  select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>% 
  gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z) 
# funcao ajuste dos modelos
f_md <- function(dados){
  var_dispersao <- unique(dados$var_dispersao)
  if(var_dispersao == "k"){
    l_md <- vector("list",2)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"))
    dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
  }else{
    dados$value_dispersao <- as.numeric(dados$value_dispersao)
    l_md <- vector("list",3)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"),
                     paste0(var_dispersao," (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         ( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    
  }
  return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio) <- NULL
l_nRef__modeloCheio <- do.call(c,l_nRef__modeloCheio)
#
f_warning <- function(glmm_object){
  # glmm_object <- l_nRef__modeloCheio[[4]]
  v_message <- glmm_object@optinfo$conv$lme4$messages %>% 
    as.character()
  if(length(v_message)==0){v_message <- "OK"}
  if(length(v_message)>1){v_message <- v_message[1] }
  return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio,f_warning,.id="glmer") %>% 
  rename(warning_message = V1)

Table 2.1 Modelos Cheios estimados e avisos de convergência

glmer warning_message
d_Lplot.z 1|Site Model failed to converge with max|grad| = 0.00343548 (tol = 0.002, component 1)
d_Lplot.z MN|Site Model failed to converge with max|grad| = 0.00836802 (tol = 0.002, component 1)
d_Lplot.z (d_Lplot.z+d_Lplot.z^2)*MN|Site OK
d.z 1|Site Model failed to converge with max|grad| = 40.1827 (tol = 0.002, component 1)
d.z MN|Site Model failed to converge with max|grad| = 0.00746139 (tol = 0.002, component 1)
d.z (d.z+d.z^2)*MN|Site OK
k 1|Site OK
k MN|Site OK

allFit

i <- 1
names(l_nRef__modeloCheio[v_glmerUpdate])[i]
#
## update1:
md_allFit <- allFit(l_nRef__modeloCheio[v_glmerUpdate][[i]],maxfun = 1e5, parallel = 'multicore', ncpus = 3)

1) d_Lplot 1|Site

7 optimizer(s) failed

2) d_Lplot MN|Site

7 optimizer(s) failed

3) d 1|Site

7 optimizer(s) failed

4) d MN|Site

7 optimizer(s) failed

Modelo cheio 2:

# dados
df_md <- df_resultados %>% 
  select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>% 
  gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z) 
# funcao ajuste dos modelos
f_md <- function(dados){
  var_dispersao <- unique(dados$var_dispersao)
  if(var_dispersao == "k"){
    l_md <- vector("list",2)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"))
    dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
  }else{
    dados$value_dispersao <- as.numeric(dados$value_dispersao)
    l_md <- vector("list",4)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"),
                     paste0(var_dispersao," ",var_dispersao,"*MN|Site"),
                     paste0(var_dispersao,"^2 (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (value_dispersao * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         ( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    
  }
  return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio2 <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio2) <- NULL
l_nRef__modeloCheio2 <- do.call(c,l_nRef__modeloCheio2)
#
f_warning <- function(glmm_object){
  v_message <- glmm_object@optinfo$conv$lme4$messages %>% 
    as.character()
  if(length(v_message)==0){v_message <- "OK"}
  if(length(v_message)>1){v_message <- v_message[1] }
  return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio2,f_warning,.id="glmer") %>% 
  rename(warning_message = V1)

Table 2.2 Modelos Cheios estimados e avisos de convergência

glmer warning_message
d_Lplot.z 1|Site Model failed to converge with max|grad| = 0.00284637 (tol = 0.002, component 1)
d_Lplot.z MN|Site OK
d_Lplot.z d_Lplot.z*MN|Site OK
d_Lplot.z^2 (d_Lplot.z+d_Lplot.z^2)*MN|Site OK
d.z 1|Site Model failed to converge with max|grad| = 0.00369235 (tol = 0.002, component 1)
d.z MN|Site OK
d.z d.z*MN|Site OK
d.z^2 (d.z+d.z^2)*MN|Site OK
k 1|Site OK
k MN|Site OK

allFit

d_Lplot.z 1|Site
7 optimizer(s) failed

d.z 1|Site
7 optimizer(s) failed

Comparação de Modelos Cheios:

Tabela 2.3 Comparação baseada em AICc dos modelos cheios

GLMM dAICc df weight
d^2 (d+d^2)*MN|Site 0.0 39 1
d_Lplot^2 (d_Lplot+d_Lplot^2)*MN|Site 274.6 39 <0.001
d d*MN|Site 38427.3 22 <0.001
d_Lplot d_Lplot*MN|Site 38586.5 22 <0.001
k MN|Site 74141.2 123 <0.001
d MN|Site 95710.0 15 <0.001
d_Lplot MN|Site 96279.6 15 <0.001
k 1|Site 106158.0 121 <0.001

Diagnostico do modelo cheio mais plausível 1

Figura 1.2 Resíduos Quantílicos do modelo cheio plausível

Diagnostico do modelo cheio mais plausível 2